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Abstract. We calculate reduced moments of the matter density fluctuations, 
up to order q — 5, from counts in cells produced by Particle-Mesh numerical 
simulations with scale-free Gaussian initial conditions. We use power-law spec- 
tra P{k) oc with indices n = —3, —2, —1, 0, 1. Due to the supposed absence 
of characteristic times or scales in our models, all quantities are expected to de- 
pend on a single scaling variable. For each model, the moments at all times can 
be expressed in terms of the variance alone. We look for agreement with 
the hierarchical scaling ansatz, according to which oc ^2'^ • -^or n < —2 
models we find strong deviations from the hierarchy, which are mostly due to 
the presence of boundary problems in the simulations. A small, residual sig- 
nal of deviation from the hierarchical scaling is however also found in n > — 1 
models. For the first time, due to our large dynamic range and careful checks 
of scaling and shot-noise effects, we are able to detect evolution away from the 
perturbation theory result. 

Subject headings: Galaxies: formation, clustering - large-scale structure of the 
Universe. 
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1. Introduction 

A fundamental problem in the analysis of the matter distribution in the 
universe is to choose simple statistical tools able to provide the most compact 
information on both the initial conditions and the subsequent evolution of den- 
sity fluctuations. In this sense, the study of count probabilities has proved a 
useful statistical technique. They can be used to follow the action of gravity 
during the mildly nonlinear as well as fully nonlinear phases of structure forma- 
tion. The count probability approach, dating back to the early work of Hubble 
(1934), has been recently applied to a number of galaxy samples: Efstathiou 
et al. (1990) calculated the variance of IRAS-selected galaxies in the QDOT 
sample for roughly cubical cells of various sizes, while Loveday et al. (1992) 
performed the same analysis in the Stromlo-APM redshift survey; Saunders 
et al. (1991) computed the skewness of density fluctuations, after smoothing 
the QDOT galaxy distribution by a Gaussian filter. A statistical analysis of 
the CfA and SSRS optical galaxy samples in terms of moments of counts in 
cells has been recently performed by Gaztanaga (1992; see also Gaztanaga & 
Yokoyama 1993). A more recent analysis of this type, up to the fifth connected 
moment, has been performed by Bouchet et al. (1993) on the 1.2 Jy IRAS 
Galaxy Redshift Survey (see also Bouchet, Davis, Sz Strauss 1992). Compared 
to connected correlation functions of order g, ^q(xi, . . . ,Xg), reduced moments 
(or cumulants) of the same order, of the fractional density fiuctuation en- 
hance the signal-to-noise ratio, though at the expense of reducing the amount 
of geometrical information. One has the following connection between the above 
quantities 

lq{R) = j (fxi...d'xqWR{^l)...WR{^q)^q{^l,...,^q), (l) 

where Wr{^) defines a suitable filter over a volume of size R. These relations 
allow one to connect the results on moments of galaxy counts in cells with the 
large amount of available data on galaxy correlation functions. Actually, mo- 
ments of counts in cells can be related to the only after shot-noise subtraction 
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(see the following Section). Early observations of higher order (i.e. q > 2) cor- 
relation functions established the validity of the so-called hierarchical scaling 
ansatz according to which correlations of order q can be expressed as suitable 
sums of products oi q — 1 two-point functions (Groth & Peebles 1977; Fry & 
Peebles 1978; Sharp, Bonometto, Sz Lucchin 1984). If the two-point function 
scales with distance as a power-law, ^(r) oc r~'^, or if the filtering radius is larger 
than the typical correlation length, a related hierarchy holds for the reduced 
moments, 

e,(i?) = %'"'(i?), 9>2, (2) 

with constant coefficients Sq. The hierarchical scaling of Eq.(2) can be given a 

theoretical justification in two different regimes. Starting from Gaussian density 

fluctuations, perturbation theory shows that the action of gravity, already in 

the mildly nonlinear regime, predicts the above hierarchical structure (Peebles 

1980; Pry 1984b; Goroff et al. 1986; Bernardeau 1992). The hierarchical scaling 

also represents a self-consistent solution of the BBGKY equations in the fully 

relaxed, highly nonlinear regime (Davis & Peebles 1977; Pry 1984a; Hamilton 

1988). The validity of this ansatz has been successfully tested in numerical 

simulations by a number of authors: Coles & Frenk (1991); Bouchet & Hernquist 

(1992, hereafter BH); Weinberg & Cole (1992, hereafter WC); Lahav et al. 

(1993, hereafter LIIS); Pry, Melott, & Shandarin (1993). Most of these works, 

however, deal with the skewness, ^3, vs. the variance, cr^ = relation (see 

also Silk & Juszkiewicz 1991). Indeed, much recent work has been devoted 

— — 2 

to a detailed analysis of the skewness ratio S3 = ^3{R)/^2 (R)- particular, 
Juszkiewicz & Bouchet (1992) and Juszkiewicz, Bouchet, & Colombi (1993) 
used a second order perturbation theory in Eulerian coordinates to compute 
the dependence of S3 on the type of window function as well as on the spectral 
index n for Gaussian density fluctuations with scale-free spectra P{k) oc k^. 
Bouchet et al. (1992) used second order perturbation expansion in Lagrangian 
space to evaluate the dependence of S3 on the density parameter Q. 
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The relevance of the primordial skewness of density fluctuations in de- 
termining both the dynamical evolution and the present texture of the matter 
distribution has been discussed by Moscardini et al. (1991), Messina et al. 
(1992) and WC. The quantity S3 has been calculated by Coles et al. (1993) for 
the mass and galaxy distribution in A'"-body simulations of skewed Cold Dark 
Matter (CDM) models and shown that this quantity can be used as a pow- 
erful test to discriminate among various statistical distributions of primordial 
fluctuations. 

From the theoretical point of view, even if we accept that the mass dis- 
tribution follows the hierarchical law on scales affected by nonlinear evolution, 
one still has to ask whether such a scaling is stable against the nonlinear (and 
possibly non-local) biasing process leading to the galaxy distribution. This is- 
sue has been partially solved by Fry & Gaztanaga (1993), who showed that the 
hierarchical scaling is indeed preserved by a rather general type of bias in the 
limit of small fluctuations. Conversely, it might be that the observed scaling of 
higher order correlation functions is entirely due to the bias mechanism, instead 
of reflecting the true statistical properties of the underlying matter distribu- 
tion. Actually, Politzer & Wise (1984), for the Gaussian case, and Matarrese, 
Lucchin, & Bonometto (1986), for the non-Gaussian one, argued that Eq.(2) is 
also implied on intermediate and large scales, if the biasing mechanism requires 
a high-density threshold, while further terms deriving from the Kirkwood ex- 
pansion (e.g. Peebles 1980) also appear on small scales. On the other hand, 
biasing with moderate threshold may lead to the hierarchical form (Melott & 
Fry 1986). Any possible contamination of the hierarchical scaling law due to 
redshift-space distortions has been shown to be negligible by a number of au- 
thors (Bouchet et al. 1992; LIIS; Coles et al. 1993; Juszkiewicz, Bouchet, & 
Colombi 1993). 

In this paper we test the validity of the hierarchical scaling law of Eq.(2) 
for the third (skewness), fourth (kurtosis) and flfth connected moments of the 
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density contrast in A'"-body simulations of the gravitational evolution of scale- 
free Gaussian models with spectral index n = —3, —2, —1, 0, 1. We also study 
the dependence of Sq (for q = 3, 4, 5) on the primordial spectral index n. 
Moments of counts in cells have already been computed in numerical simulations 
for some of these models. Bouchet & Hernquist (1992), besides considering 
CDM and Hot Dark Matter models, run a tree code with white-noise (i.e. n = 
0) initial conditions. In the frame of a comparison of Particle-Mesh simulations 
with both Gaussian and non-Gaussian initial conditions, WC calculate for 
n = —2, —1, initial spectra. Finally, LIIS analyze counts in cells in tree code 
simulations for various models, including n = —1, 0, 1 scale-free models. Our 
analysis, besides considering a wider ensemble of power-law models, covers a 
much larger dynamic range both in time and resolution. In Table I we show the 
range in the rms fluctuation a and spectra studied for our simulation and those 
cited in this paragraph. It appears that ours is the only study to date with 
suflBcient dynamic range and control of boundary conditions to reliably detect 
evolution away from the perturbation theory result. We also have looked at a 
variety of pure power-law models, so that such dependence can be detected. 

The plan of the paper is as follows. In Section 2 we give the theoretical 
background. Section 3 presents the numerical simulations and the cautions used 
in the following analysis, while Section 4 discusses our results on the analysis 
of the moments of counts in cells. Conclusions are summarized in Section 5. 

2. Moments of Counts in Cells 

The quantities ^q{R) defined above represent the reduced moments of the 
density fluctuation (5i?(x) = qii{-x)/'q — 1 (here qh is the density smoothed over 
the scale R and 'g its average) ; these can be derived from the moment generating 
function M{s) = J dgnVign) exp{isgii/'g)^ where V{gR) gives the probability 
density of the continuous variable gR{'x). One has 




q=2 



(3) 
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Given the probability density or the moment generating function, one can easily 
generate discrete count-probabilities Vm (e.g. Peebles 1980). In fact, the count 
Vm can be understood as the probability that, in a realization of the stochastic 
process qr, m objects are found in a randomly placed cell of volume V = Nv, 
where v is the specific volume 1 /'g and N is the expected number of objects in 
that volume. Consider then gji as giving the mean density of an ensemble of 
local Poisson count distributions with mean Nvgn, whose moment generating 
function is cxp[A^f ^ij(e*^ — 1)]. Note that the Poisson model does not necessarily 
provide a good representation of discreteness effects (see e.g. Coles & Frenk 
1991; Borgani et al. 1993); in particular, it is likely to fail when the expected 
number of objects N in the cell volume is smaller than unity, i.e., when V <^ v. 
Averaging over the gji ensemble produces the moment generating function for 
the discrete process as 



Mdisis) = / dgR V{gR) exp[NvgR{e'' - 1)] = M[-iie'' - 1)N] : (4) 



the moment generating function of discrete counts, J^dis{s), is obtained from 
the continuous one, M.{s), by the replacement is (e** — 1)A'". By inverse 
Fourier transforming Mdisis) one gets VdisiQu) = Z)m=o KQR-'^Q)'Pm, where 
5 is the Dirac delta-function and the counts Vm are defined through a Poisson 
transform of V{gR), 



One can then define the normalized central moments of the counts in cells 
of volume V as Hq = {{{m — m)/m)'^), where 



and fn = ^^^^ wP^n = N. For small A^, at fixed moments of gR/'g, the 
count distribution is shot-noise dominated and Addis (s) reduces to the moment 
generating function of a Poisson process with mean N: if the cell volume V is 
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too small, statistical fluctuations dominate the realization and one is unable to 
get any faithful statistical information on Vigu) from the counts. This can be 
seen by writing 'P{qr) in Eq.(4) through its Fourier transform and expanding 
in reduced moments ^q{R), 



q=2 

which for small N and fixed $,q{R) yields 

Mdisis) ^ j dy exp[y(e^^ - 1)] 5{y - iV) = exp[iV(e^^ - 1)]. 



, (7) 



(8) 



When dealing with the — > limit above, one should, however, consider the 
volume dependence of the connected moments. Actually, if the hierarchical 
scaling of Eq.(2) holds and oc N~'^/^ as A?" — > 0, shot-noise dominates for 
small filtering scales, provided that the effective spectral index ne// = 7 — 3 is 
smaller than zero. Note that on small scales n^ff 7^ n. 

Conversely, shot-noise may even dominate for large cell sizes, N ^ 00, 
where 7 — > n + 3, if the spectral index n is larger than zero. Except for these 
cases, one generally expects that the discrete counts NVm reduce to the original 
continuous distribution of qh, for large N and m and fixed m/N, as a property 
of the Poisson transform. In fact, using the asymptotic representation of the 
Poisson counts. 



[NvQRr -Nvsn 1 

e ~ . = exp 

ml y/2'KNvgR 



(m - NvQ^f 
2NvgR 



(9) 



in the integrand of Eq.(5), and taking the limit for AT —> 00 one gets NVm ~ 
Vign), with qr = gm/N. 

We can conclude that the optimal range of cell sizes R = yV3 depends 
on the spectral slope n: in a numerical simulation such as ours, with mean 
interparticle distance £, one should require R > Max[^, £ for every 

n; for n > 0, however, one should also require R < £ c7~^/^(i?). 



From Eqs.(3), (4) and (6) above one can explicitly find the required re- 
lations among the moments of the counts and the reduced moments of the 
continuous variable 6r. Inverting these relations, one finally obtains the cumu- 
lants as a function of the moments of counts in cells /in, up to order q. We 
have, in particular, 

l2 = /'2-^, (10) 
{3 = W - 3^ + ^, (11) 

f4 = /'4-6^-3,.i+ll^-A (12) 
= - lof - (10^. - §)^3 + 30| - 50^ + ^. (13) 

— — 2 — — 3 

In what follows we shall also consider the ratios ^'3 = C3/C2 ) 'S'4 = C4/C2 

— — 4 

and S5 = C5/C2 ; order to test whether the hierarchical scaling relations 
apply, i.e., whether these ratios are scale-independent, i.e. independent of the 
variance. 

Theoretical predictions for the value of S3 have been obtained using second 
order perturbation theory. In order to get a consistent prediction for S4, and S5 
one respectively needs third and fourth order perturbative calculations. Goroff 
et al. (1986) have computed these ratios for initially Gaussian perturbations in 
standard CDM, by a clever summation of tree diagrams. Filtering the density 
field by a Gaussian window they obtain the values Ss ^ 3, S4 !v 16 and ^'5 ~ 100 
on large scales, where the spectral slope tends to the Zel'dovich value n = 1. 
For scale-free Gaussian initial perturbations in an Einstein-de Sitter model, 
Juszkiewicz, Bouchet, & Colombi (1993) find the relation 

'5?^ = y-(n + 3), -3<n<l, (14) 

using a spherical top-hat filter, while for n = 1 the perturbative prediction 
formally diverges. In the latter case, taking into account that numerical sim- 
ulations cannot reproduce the initial spectrum above the Nyquist frequency, 
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they find S^^' {n = 1) = 1.9. We shall compare these perturbative estimates 
with our numerical results. The use of a sharp cubic filter instead of a spherical 
one is not expected to introduce big changes in the 5*3 vs. n relation. Actu- 
ally, we have numerically verified that the two filters give essentially equivalent 
results, provided all quantities are compared at equal smoothing volumes (see 
also, LIIS). 

Perturbation theory implies an expansion of a series. As is well known 
in basic physics, the series contains higher and higher order powers of the per- 
turbed quantity. It can only be expected to converge to the correct result if this 
quantity is small. Indeed, perturbation theory is going to fail as a gets of order 
unity, simply because the gravitational field becomes arbitrarily large around 
regions of orbit mixing. Thus perturbation theory ought to give better results 
for small a, and we can use that to compare with our procedures. The value 
of numerical simulations such as these is that we can investigate the nonlinear 
regime. 

We believe that, for the first time, due to our large dynamic range and 
careful checks of scaling, we are able to detect the evolution away from the 
perturbation theory result. 

3. Numerical Simulations 

The simulations studied herein are numerical models of the clustering of 
coUisionless matter in an expanding universe. We wish to investigate the above 
scale-invariant behavior in the case of Gaussian initial conditions. In order 
to implement this, we use an O = 1 universe, as to choose otherwise would 
introduce a preferred scale or time. We use pure power-law initial perturbation 
spectra P{k) oc /c" with n = —3, —2, —1, 0, 1. 

The simulations are done with a Particle-Mesh (PM) code (Hockney & 
Eastwood 1981) with 128^ particles in as many cells. In this paper we use 
10 runs (two of each spectral index) out of an ensemble of 50 generated for 
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other systematic studies (Melott & Shandarin 1993). The PM code used in 
these studies is highly optimized, using a staggered mesh scheme, and has twice 
the dynamical resolution of any other PM code with which it has been com- 
pared [Melott 1986; Melott, Weinberg, & Gott 1988 (hereafter MWG); Wein- 
berg 1993a,b]. Thus the studies shown here are roughly equivalent in dynamic 
range to usual PM runs with 128^ particles on a 256^ mesh except that we have 
less shot-noise and coUisionality. Having a relatively large number of particles 
has the advantage of good mass resolution and the ability to impress initial per- 
turbations right up to the particle Nyquist frequency k^y = 64. Other methods 
such as P^M and tree codes have not yet been able to run with 128^ particles, 
which is relatively easy with PM, as can be seen by our large ensemble of such 
runs. More details about the particular simulations used here can be seen in 
Melott & Shandarin (1993). 

Having stressed some advantages of our simulations, we would now like to 
discuss some of the precautions needed in using them, particularly for studies of 
scale-free processes. Resolution is one problem. Resolution has been stressed as 
an advantage of codes in which short-range forces are calculated accurately be- 
tween point masses. In reality there are a number of different but related kinds 
of resolution: (a) mass resolution, essentially the reciprocal of the number of par- 
ticles; (b) force resolution, essentially how accurately the force law tracks 1/r^ 
at small separations; (c) spectral resolution, equivalent to the minimum of the 
number of particles or Fourier analysis cells per unit length, whichever is smaller; 
(d) minimum scale on which two-body relaxation becomes important. PM 
methods are superior for all except (b), in which other methods work better. 
We first consider this limitation. 

The growth rate of various modes in linear theory was studied in MWG 
for this PM code. The growth rates for PM codes are usually described as 
being unacceptable for k > 0.25 /cjvy, or equivalently for wavelength less than 
8 cells. MWG found an equivalent performance for 0.5 kNy to that found in 
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usual PM codes at 0.25 k^y This results from the staggered mesh scheme; it 
has since been confirmed in other comparisons (Weinberg 1993a,b) and the PM 
code used here has been slightly improved since that corresponding to a further 
30% resolution increase, i.e. giving results at A = 3 cells equivalent to usual 
PM codes at A = 8 cells. 

We still must take account of limited force resolution based on our grid 
scheme. The advantage of the methods used here is that we can test for the 
adequacy of our precautions by observing the results at various stages. In 
pure power-law models, properties of the distribution should depend only on 
a = {Sg/g)rms, assuming the use of the identical initial power-law smoothing 
windows. But a increases with time, and the agreement on difi'erent scales at 
different times with the same u is a strong consistency test. A similar strategy 
has already been used to find previously unknown effects from the absence of 
waves larger than the simulation volume (Kauffmann & Melott 1992; Gramann 
1992; Melott & Shandarin 1993). 

We first describe the stages of our simulations and then the restrictions 
we applied. Our simulations were started by using the Zel'dovich (1970) ap- 
proximation, as first utilized by Klypin & Shandarin (1983). It is well-known 
that this approximation is inaccurate beyond the time of nonlinearity [although 
better than other approximations studied, with appropriate filtering; see Coles, 
Melott, Sz Shandarin (1993)] so the initial amplitude is well below unity at the 
Nyquist frequency. 

The simulations were stopped at kne = kNy, 0.5 kNy, 0.25 A;jvy, . . . , 2 A;/ 
where kf = 2tt/L is the fundamental mode of the box; kni is defined by 



In this study we have an available range of kni from 2 /c/ to 64 /c/. Everything in 
the simulations should scale as l/kn£, in the absence of boundary or resolution 




(15) 
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problems. We make the following restrictions on what scales will be studied. We 
will show results based on counts in cells, in boxes of various sizes, at various 
stages of nonlinearity. All results for a given spectral index will be plotted 
together as a function of c — {5q/ Q)rms 

and reveal any problems. 

(a) The stage with kni = kNy will not be used to study nonlinear effects since 
the code is known not to perform well at this frequency. We will begin 
with the stage kne = 0.5 /cTVy = 32 /c/ in this study, except that we use 
the earlier stage to help establish the linear limit. 

(b) At each stage we will not present results for one pixel of density, since this 
depends primarily on kNy The smallest will be a cube of 8 such cells. We 
expect this will be acceptable since our code performs well at wavelength 
of 4 cells, and collapse of A = 4 perturbations will initially give rise to 
condensations of diameter 2. We will confirm this later. 

(c) Kauffmann & Melott (1992) found that for voids of size greater than 
size 0.25L self-similarity was broken in a model equivalent to our index 
n = -1; see also Gramann (1992) and Melott & Shandarin (1993). We 
therefore restrict ourselves to cubes of size L/8 or smaller. Combined with 
restriction (b), this leaves us with cells for scale from L/64 to L/8. We do 
not expect this to work with n < —2, and believe that for these values all 
so-called "N-body" simulations are at best crude. It would not be entirely 
a joke to say that in this case a model could never be big enough to be a 
fair sample of itself. Fortunately, it appears that the power spectrum of 
the universe turns over to n > on large scales. 

(d) We use counts in cells, rather than the usual Cloud-in-Cell method (Hock- 
ney & Eastwood 1981) to bin densities. This procedure increases shot- 
noise as compared to that present when the PM code calculates the grav- 
itational potential. For this reason we do not use cell sizes with a < 0.1. 
In practice, in our simulations this restriction eliminates cells where the 
shot-noise power is comparable to the fluctuation power impressed on the 
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simulation, which makes the subtraction doubtful. Following our discus- 
sion in Section 2, this only affects the largest cells at early times in our 
models with n > 0. 

Taken together, these restrictions can eliminate problems and give us much 
more usable dynamic range than has been possible before in such a study. 
We verify that it works remarkably well by examining the agreement between 
various stages at the same cr. It suggests that scale dependence found in previous 
studies was probably a result of boundary problems. For n < — 2 we find very 
poor agreement, as expected. 

As an illustration of the different appearance of the models, Figure 1 
shows, for the five spectra, grey-scale plots of L/64 thickness slices at the stage 
corresponding to knt = 8. 

4. Analysis and Results 

In our scale-free simulations we expect that every quantity depends on a 
single scaling variable, namely the variance ^2- Therefore we can plot the results 
of counts in cells at all times together as a function of the variance. Figure 2 
shows the skewness ^3 vs. the variance, for all models, i.e. for all values of n. 
The points shown are the average of the two runs; the corresponding dispersion 
is always small and will not be shown here. The dashed lines are the second 
order perturbative predictions [obtained from Eq.(14)], for the same value of 
n; these computations derive from using a spherical top-hat window, but we 
checked that changing from spherical to cubic filter produces essentially the 
same results. The solid line shows the two-parameter linear fit 

log^^^A, + B,\og^^ (16) 

(for q = 3); the corresponding coefficients and related errors are reported in 
Table II. Figures 3 and 4, respectively, show the kurtosis ^4 and the fifth moment 
^5 vs. for all values of the initial spectral index. Perturbative predictions 
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are not available in these cases. The solid lines represent the results of linear 
fits from Eq.(16); best-fit coefficients and errors are reported in Table III and 
IV, for the fourth and fifth moment, respectively. Note that the n < —2 results 
are clearly not reliable, being highly affected by the finite box-size (see Fry, 
Melott, Sz Shandarin 1993). In all other cases (n = —1, 0, 1), the scaling of 
{q = 3, 4, 5) with the variance is quite close to the hierarchical form, Bq = q—1. 
Nevertheless, we detect a residual deviation from this scaling, at more than 
three standard deviations for both the skewness and the kurtosis, while for the 
fifth moment this result appears at only one standard deviation. We consider 
rather unlikely that such a deviation is produced by boundary effects. Due 
to the large dynamic range of our simulations and to the careful treatment 
of resolution and shot-noise problems, we are led to trust these deviations, 
even though our results are closer to the hierarchy than previous works (BH 
and LIIS), whose simulations are more affected than ours by finite sample and 
resolution effects. For instance, for initial white-noise, BH find ^3 = 2.10±0.01, 
B4 = 3.26 ± 0.02 and B5 = AAA ± 0.04, while LIIS obtain S3 = 2.08 ± 0.01 and 
-B4 = 4.16 ± 0.03. Note however, that, according to our previous discussion on 
shot-noise, the highest a points appearing in BH, which come from cells with 
size smaller than the mean interparticle distance, are likely to be dominated by 
discreteness effects, which decreased their statistical reliability. 

In order to better display possible deviations from scale-invariance we also 
plot the coefficients Sq as a function of cr. In Figures 5, 6 and 7 we show S'3, 
S4 and S5, for all values of n. The meaning of the solid lines and dashed ones 
(when present) is as before. Note that the trend of Sq with n is the same for all 
q: Sq decreases with increasing n. Although this qualitative trend is the same 
predicted by perturbation theory for q — 3, the value of the coefficients is only 
in partial agreement with it. In Tables II - IV, we report the coefficients of the 
one-parameter fits of Sq, obtained at fixed hierarchical slope; dotted lines in 
Figures 5-7 represent these best-fit coefficients. Note that the values of S^ are 
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generally different from as given by Eq.(14); however, since perturbation 
theory is at most consistent with mildly nonlinear evolution, we also estimated 
the skewness ratio by fitting only points with a < 0.7: this is reported as in 
Table II (we similarly define ^4 and ^'5 in Tables III and IV). The agreement 
with perturbation theory is indeed improved, but we still get slightly higher 
values: = 4.5, 3.4, 3.1, 2.0, 1.7, instead of S^^^ = 4.9, 3.9, 2.9, 1.9, 1.9, 
for n = —3, —2, —1, 0, 1, respectively. In deriving the values of 5*3 above, we 
also used an earlier stage of the simulations corresponding to k^i — 64 kf. For 
their simulations with n = 0, BH argue that ^3 should tend to about 1.8 at 
large scales, very close to the perturbative prediction; their result is obtained 
from large cell sizes, corresponding to about one quarter of the computational 
box size, where sample effects make the resulting data less reliable. 

Let us stress that agreement with second order perturbation theory is not 
a good test of the quality of A'"-body data, except in the linear regime. However, 
the behavior of such codes has already been widely tested in this regime. In 
the nonlinear regime, one should use A'"-body data to test the reliability of 
perturbation expansion results, such as the second order estimate of Ss reported 
in Eq.(14). In this sense one can consider a success that the qualitative trend of 
S3 with n is correctly predicted by perturbation theory! It would be interesting 
to have similar predictions for higher order moments, to compare with our 
numerical results. 

5. Conclusions 

We have studied the properties of higher order moments of the matter 
distribution generated by gravity in the nonlinear regime. This was done by 
analyzing numerical simulations of the evolution of initially Gaussian perturba- 
tions with scale-free power-spectra (with spectral index n = —3, —2 , —1, 0, 1) 
at many evolution stages and smoothing scales. Our results for n > — 1 models 
indicate that these moments are fairly close to the hierarchical scaling ansatz, 
according to which connected moments of order q, are proportional to the 
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q — 1 power of the variance However, we detect a residual dependence, 
above the statistical noise, of the coeflBcients Sq = iq/ ^2 ^ on the variance, 
i.e. on scale. In order to detect such a signal from the data we had to properly 
account for shot-noise, finite resolution effects and boundary problems. For 
models with n < —2, where the amount of large-scale power does not allow a 
fair representation of low frequency modes of the density field (due to the finite 
box size), a stronger dependence of Sq on scale is found. 

The skewness ratio ^'3 is found to decrease with increasing spectral slope, 
n, i.e. with decreasing large-scale power, as correctly predicted by perturbative 
calculations, although our values for are only consistent with the theory if 
the small a limit of these quantities is considered. The same qualitative trend 
with n is seen for the higher order coefficients 5*4 and 6*5. 

Altogether, these results indicate that: 1) the dynamical effect of gravity 
is such as to generate non-Gaussian signatures on a Gaussian initial density 
field, already in the earliest stages of evolution and/or on large scales; 2) the 
hierarchical ansatz provides only an approximate description of the behavior of 
higher order moments of the density fiuctuation field with the scale. 
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Figure captions 



Figure 1. Grey-scale plots of the projected density in slices of thickness L/64 
at stage kn(, = 8k f in our simulations. Regions below the mean density are 
white; regions of density p > 10 are black; a grey scale is used in between, (a) 
n = -3; (b) n = -2; (c) n = -1; (d) n = 0; (e) n = 1. 

Figure 2. The skewness vs. the variance ^2 for the different models. The 
solid lines represent the best-fit obtained from Eq.(16), while the dashed lines 
are the perturbative prediction. Different symbols refer to different stages of 
the simulations: k^e = 32 kf, filled triangles; kni = 16 kf, asterisks; k^e = 8 kf, 
open triangles; knt = 4 kf, open squares; kne = 2 kf, open circles. Note that 
circles are absent for n = —3; they were left out due to lack of any benefit in 
including them. 

Figure 3. The kurtosis ^4 vs. the variance ^2 foi" the different models. The 
solid lines represent the best-fit obtained from Eq.(16). Different symbols refer 
to different stages of the simulations as in Figure 2. 

Figure 4. The fifth connected moment ^5 vs. the variance ^2 for the different 
models. The solid lines represent the best-fit obtained from Eq.(16). Different 
symbols refer to different stages of the simulations as in Figure 2. 

— — 2 

Figure 5. The skewness coefficient S3 = vs. the rms fluctuation a 

for the different models. The solid lines represent the best-fit obtained from 

Eq.(16), the dashed lines are the perturbative prediction, finally the dotted lines 

are the result of a best-fit forced to the hierarchical slope. Different symbols 

refer to different stages of the simulations as in Figure 2. 

_ 3 

Figure 6. The kurtosis coefficient S4 = ^4/^2 vs. the rms fluctuation a for the 
different models. The solid lines represent the best-fit obtained from Eq.(16), 
while the dotted lines are the result of a best-fit forced to the hierarchical slope. 
Different symbols refer to different stages of the simulations as in Figure 2. 
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Figure 7. The fifth moment coefficient S5 = ^^e rms fluctuation 

a for the different models. The solid lines represent the best-fit obtained from 
Eq.(16), while the dotted lines are the result of a best-fit forced to the hierar- 
chical slope. Different symbols refer to different stages of the simulations as in 
Figure 2. 
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Table I 

Comparison of Dynamic Range in Recent Similar Studies 

Ratio of amax to amin 



BH* WC** LIIS*** this study 
-3 ~ 10 

-2 1.86 ~ 20 

-1 1.86 ~ 8 ~ 64 

~ 100 1.86 ~ 15 ~ 100 

+1 ^ 30 ~ 100 



* estimated from BH Figure 8. 

* estimated from WC Table 2. WC did not attempt to study the evolution 
of these moments over a wide dynamic range. 

* numbers estimated from LIIS Figure 2. Discrepancies up to a factor two 
exist between results at different moments (see LIIS Table I). 



Table II 
Third moment coefRcients 





^3 


S3 


-53 


St 


q(p) 
'-'3 


-3 


0.74 ±0.01 


2.28 ±0.02 


7.0 


4.5 


4.9 


-2 


0.57 ±0.01 


2.07 ±0.01 


4.0 


3.4 


3.9 


-1 


0.51 ±0.01 


2.03 ±0.01 


3.3 


3.1 


2.9 





0.38 ±0.01 


2.04 ±0.01 


2.5 


2.0 


1.9 


+1 


0.30 ±0.01 


2.04 ±0.01 


2.0 


1.7 


1.9# 



see discussion in the text. 
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Table III 
Fourth moment coefficients 



n 


-3 


1.66 ±0.02 


3 


II 


-2 


1.33 ±0.02 


3 


II 


-1 


1.19 ±0.02 


3 


II 





0.96 ±0.02 


3 


II = 


+1 


0.73 ±0.02 


3 



B4 S4 S'4 

59 ± 0.05 97.4 29.5 

15 ±0.03 28.1 19.1 

07 ±0.02 18.1 16.0 
04 ±0.01 10.0 9.4 

08 ±0.02 6.2 5.9 



Table IV 
Fifth moment coefficients 







A, 


55 


-55 


St 


n 


-3 


2.65 ±0.04 


4.92 ± 0.08 


2091 


218 


II 


-2 


2.15 ±0.04 


4.24 ±0.05 


319 


161 


II 


-1 


2.00 ± 0.03 


4.05 ±0.03 


147 


135 


II 





1.64 ±0.03 


4.02 ±0.03 


63 


86 


II = 


±1 


1.31 ±0.03 


4.02 ±0.03 


30 


48 
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